#!/bin/bash -e
# based on ug.sh, filter is more easier.
function info() {
echo Usage: `basename $0` '[-i intervals.bed] [-G(somatic)1(single_cell)m(mixture)s(shallow_coverage_between_30_100)] [-d depth_design] in.bam'
exit 65
}

while getopts  ":i:p:G1smd:" opts
do
        case  $opts  in
        i) intervals=$OPTARG;;
		p) out_prefix=$OPTARG;;
		G) somatic=true;;
		1) single_cell=true;;
		m) cells_mixture=true;;
		s) shallow_seq=true;;
        d) depth_design=$OPTARG;;
		*) info;;
        esac
done
shift $(($OPTIND - 1))

if test -z "$1"; then info; fi


. $var

nt=3
nct=2
mem=8

depth_design=${depth_design:=30}

# --min_base_quality_score / -mbq
# this argument is ignored in indel calling. In indel calling, low-quality ends of reads are clipped off (with fixed threshold of Q20).

# --min_indel_count_for_genotyping / -minIndelCnt
# Minimum number of consensus indels required to trigger genotyping run
# A candidate indel is genotyped (and potentially called) if there are this number of reads with a consensus indel at a site. Decreasing this value will increase sensitivity but at the cost of larger calling time and a larger number of false positives.

# --min_indel_fraction_per_sample / -minIndelFrac
# Minimum fraction of all reads at a locus that must contain an indel (of any allele) for that sample to contribute to the indel count for alleles
# Complementary argument to minIndelCnt. Only samples with at least this fraction of indel-containing reads will contribute to counting and overcoming the threshold minIndelCnt. This parameter ensures that in deep data you don't end up summing lots of super rare errors up to overcome the 5 read default threshold. Should work equally well for low-coverage and high-coverage samples, as low coverage samples with any indel containing reads should easily over come this threshold.

# --standard_min_confidence_threshold_for_calling / -stand_call_conf
# The minimum phred-scaled confidence threshold at which variants should be called
# The minimum phred-scaled Qscore threshold to separate high confidence from low confidence calls. Only genotypes with confidence >= this threshold are emitted as called sites. A reasonable threshold is 30 for high-pass calling (this is the default).

# --standard_min_confidence_threshold_for_emitting / -stand_emit_conf
# The minimum phred-scaled confidence threshold at which variants should be emitted (and filtered with LowQual if less than the calling threshold)
# This argument allows you to emit low quality calls as filtered records.

# --pcr_error_rate / -pcr_error
# The PCR error rate to be used for computing fragment-based likelihoods
# The PCR error rate is independent of the sequencing error rate, which is necessary because we cannot necessarily distinguish between PCR errors vs. sequencing errors. The practical implication for this value is that it effectively acts as a cap on the base qualities.

# --annotateNDA / -nda
# If provided, we will annotate records with the number of alternate alleles that were discovered (but not necessarily genotyped) at a given site
# Depending on the value of the --max_alternate_alleles argument, we may genotype only a fraction of the alleles being sent on for genotyping. Using this argument instructs the genotyper to annotate (in the INFO field) the number of alternate alleles that were originally discovered at the site.

# --comp / -comp # useful for some annotation
# Comparison VCF file
# If a call overlaps with a record from the provided comp track, the INFO field will be annotated as such in the output with the track name (e.g. -comp:FOO will have 'FOO' in the INFO field). Records that are filtered in the comp track will be ignored. Note that 'dbSNP' has been special-cased (see the --dbsnp argument).
# This argument supports reference-ordered data (ROD) files in the following formats: BCF2, VCF, VCF3

# --contamination_fraction_to_filter / -contamination
# Fraction of contamination in sequencing data (for all samples) to aggressively remove
# If this fraction is greater is than zero, the caller will aggressively attempt to remove contamination through biased down-sampling of reads. Basically, it will ignore the contamination fraction of reads for each alternate allele. So if the pileup contains N total bases, then we will try to remove (N * contamination fraction) bases for each alternate allele.

# --genotyping_mode / -gt_mode
# Specifies how to determine the alternate alleles to use for genotyping
# The --genotyping_mode argument is an enumerated type (GenotypingOutputMode), which can have one of the following values:
# DISCOVERY
# The genotyper will choose the most likely alternate allele
# GENOTYPE_GIVEN_ALLELES
# Only the alleles passed by the user should be considered.

# --indelGapContinuationPenalty / -indelGCP
# Indel gap continuation penalty, as Phred-scaled probability. I.e., 30 => 10^-30/10
# byte  10

# --pcr_error_rate / -pcr_error
# The PCR error rate to be used for computing fragment-based likelihoods
# The PCR error rate is independent of the sequencing error rate, which is necessary because we cannot necessarily distinguish between PCR errors vs. sequencing errors. The practical implication for this value is that it effectively acts as a cap on the base qualities.

# --baqGapOpenPenalty / -baqGOP
# BAQ gap open penalty
# Phred-scaled gap open penalty for BAQ calculation. Although the default value is 40, a value of 30 may be better for whole genome call sets.

min_confidence_calling=15
min_confidence_emit=15
cov_max=10000 # cells_mixture needed, because min_indel_fraction=0.01, min deepth is 100*2
contamination=0.001
mbq=20 # after bqrc, mbp will decrease in general, so we use 20 instead of 30
# genotype_mode=GENOTYPE_GIVEN_ALLELES
baqGOP=40.0
pcr_error=0.0001


if test "$cells_mixture" = "true"; then
	min_indel_count=2 # 2*100*2=400, so in mixture mode, min deepth is 400
	min_indel_fraction=0.01 # .25 default, need no less than said if 2 * 100 * 2
	max_alleles=6
	pcr_error=0.0005
    if test $depth_design -lt 200; then
        min_indel_count=3
    fi
else # single cell or bulk
	min_indel_count=3 # 5*10*2=100
	min_indel_fraction=0.1 # min deepth is 100
	max_alleles=3
    if test $depth_design -lt 30; then
        min_indel_count=4
    fi
fi

if test "$single_cell" = "true"; then
	pcr_error=0.001
    min_indel_fraction=0.05
fi

if test "$shallow_seq" = "true"; then
	min_indel_count=3
fi

if test "$wg" = "true"; then
	baqGOP=30.0 # whole genome as default
fi


if test -z "$intervals"; then

echo gatk UnifiedGenotyper snp recommand nt 6 nct 4 mem 16 or "3/8/32"
java $j_mem -jar $gatk \
	-T UnifiedGenotyper \
	-R $ref_genome \
	-I $1 \
	-o $out_prefix.snp.vcf \
	-glm SNP \
	-baq CALCULATE_AS_NECESSARY \
	--baqGapOpenPenalty $baqGOP \
	--downsample_to_coverage $cov_max \
	-nt $nt \
	-nct $nct \
	-stand_call_conf $min_confidence_calling \
	-stand_emit_conf $min_confidence_emit \
	--min_base_quality_score $mbq \
	-contamination $contamination \
	--max_alternate_alleles $max_alleles \
	--annotation VariantType \
	--annotateNDA \
	--pcr_error_rate $pcr_error


echo gatk UnifiedGenotyper indel recommand nt 6 nct 4 mem 16 or "3/8/32"
if ! $(java $j_mem -jar $gatk \
	-T UnifiedGenotyper \
	-R $ref_genome \
	-I $1 \
	-o $out_prefix.indel.vcf \
	-glm INDEL \
	-baq CALCULATE_AS_NECESSARY \
	--baqGapOpenPenalty $baqGOP \
	--downsample_to_coverage $cov_max \
	-nt $nt \
	-nct $nct \
	--min_indel_fraction_per_sample $min_indel_fraction \
	--min_indel_count_for_genotyping $min_indel_count \
	-stand_call_conf $min_confidence_calling \
	-stand_emit_conf $min_confidence_emit \
	-contamination $contamination \
	--max_deletion_fraction 0.2 \
	--max_alternate_alleles $max_alleles \
	--annotation VariantType \
	--annotateNDA \
	--pcr_error_rate $pcr_error)
# 0.2*100=20 0.2*20=4
then
java $j_mem -jar $gatk \
	-T UnifiedGenotyper \
	-R $ref_genome \
	-I $1 \
	-o $out_prefix.indel.vcf \
	-glm INDEL \
	-baq CALCULATE_AS_NECESSARY \
	--baqGapOpenPenalty $baqGOP \
	--downsample_to_coverage $cov_max \
	-nt $nt \
	-nct 1 \
	--min_indel_fraction_per_sample $min_indel_fraction \
	--min_indel_count_for_genotyping $min_indel_count \
	-stand_call_conf $min_confidence_calling \
	-stand_emit_conf $min_confidence_emit \
	-contamination $contamination \
	--max_alternate_alleles $max_alleles \
	--annotation VariantType \
	--annotateNDA \
	--pcr_error_rate $pcr_error
fi

else
echo gatk UnifiedGenotyper snp recommand nt 6 nct 4 mem 16 or "3/8/32"
java $j_mem -jar $gatk \
	-T UnifiedGenotyper \
	-R $ref_genome \
	-I $1 \
	-o $out_prefix.snp.vcf \
	-glm SNP \
	-baq CALCULATE_AS_NECESSARY \
	--baqGapOpenPenalty $baqGOP \
	--downsample_to_coverage $cov_max \
	-nt $nt \
	-nct $nct \
	-stand_call_conf $min_confidence_calling \
	-stand_emit_conf $min_confidence_emit \
	--min_base_quality_score $mbq \
	-contamination $contamination \
	--max_alternate_alleles $max_alleles \
	-L ${intervals} \
	--annotation VariantType \
	--annotateNDA \
	--pcr_error_rate $pcr_error



echo gatk UnifiedGenotyper indel recommand nt 6 nct 4 mem 16 or "3/8/32"
#deal with Somehow the requested coordinate is not covered by the read?
if ! $(java $j_mem -jar $gatk \
	-T UnifiedGenotyper \
	-R $ref_genome \
	-I $1 \
	-o $out_prefix.indel.vcf \
	-glm INDEL \
	-baq CALCULATE_AS_NECESSARY \
	--baqGapOpenPenalty $baqGOP \
	--downsample_to_coverage $cov_max \
	-nt $nt \
	-nct $nct \
	--min_indel_fraction_per_sample $min_indel_fraction \
	--min_indel_count_for_genotyping $min_indel_count \
	-stand_call_conf $min_confidence_calling \
	-stand_emit_conf $min_confidence_emit \
	-contamination $contamination \
	--max_deletion_fraction 0.2 \
	--max_alternate_alleles $max_alleles \
	-L ${intervals} \
	--annotation VariantType \
	--annotateNDA \
	--pcr_error_rate $pcr_error)
then
java $j_mem -jar $gatk \
	-T UnifiedGenotyper \
	-R $ref_genome \
	-I $1 \
	-o $out_prefix.indel.vcf \
	-glm INDEL \
	-baq CALCULATE_AS_NECESSARY \
	--baqGapOpenPenalty $baqGOP \
	--downsample_to_coverage $cov_max \
	-nt $nt \
	-nct 1 \
	--min_indel_fraction_per_sample $min_indel_fraction \
	--min_indel_count_for_genotyping $min_indel_count \
	-stand_call_conf $min_confidence_calling \
	-stand_emit_conf $min_confidence_emit \
	-contamination $contamination \
	--max_deletion_fraction 0.2 \
	--max_alternate_alleles $max_alleles \
	-L ${intervals} \
	--annotation VariantType \
	--annotateNDA \
	--pcr_error_rate $pcr_error
fi

fi

. $cmd_done
